;; Naive B-spline functions
;; Copyright (C) 2023  Rocks Mazama
;;
;; This program is free software: you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation, either version 3 of the License, or
;; (at your option) any later version.
;;
;; This program is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;; GNU General Public License for more details.
;;
;; You should have received a copy of the GNU General Public License
;; along with this program.  If not, see <https://www.gnu.org/licenses/>.

(ns mazama.compgeom.b-spline)

(defn- omega [knots i p x]
  (let [tᵢ (nth knots i)
        tᵢ₊ₚ (nth knots (+ i p))
        denom (- tᵢ₊ₚ tᵢ)]
    (if (zero? denom)
      0
      (/ (- x tᵢ) denom))))

(defn- B [knots i p x]
  (if (<= p 0)
    (if (and (<= (nth knots i) x) (< x (nth knots (inc i))))
      1
      0)
    (+ (* (omega knots i p x)
          (B knots i (dec p) x))
       (* (- 1 (omega knots (inc i) p x))
          (B knots (inc i) (dec p) x)))))

(defn- naive
  "Evaluate the B-spline at `x` for degree `p`, `control-points`, and `knots`.
  `f+` is a function that adds two points together;
  `f*` is a function that multiplies a scalar and a point."
  [f+ f* control-points knots p x]
  (->> control-points
       (map-indexed #(f* (B knots %1 p x) %2))
       (reduce f+)))

(defn- lerp
  "Parametric linear interpolation from `p0` to `p1`.
  `f+` is a function that adds two points together;
  `f*` is a function that multiplies a scalar and a point."
  [f+ f* p0 p1 t]
  (f+ (f* (- 1 t) p0) (f* t p1)))

(defn- alpha [x tₖ₊₁ tᵢ]
  (if (== tₖ₊₁ tᵢ)
    0
    (/ (- x tᵢ) (- tₖ₊₁ tᵢ))))

(defn- step [f+ f* d x tₖ₊₁ tᵢ]
  (let [⍺ (map #(alpha x %1 %2) tₖ₊₁ tᵢ)]
    (map #(lerp f+ f* %1 %2 %3) d (rest d) ⍺)))

;;; https://en.wikipedia.org/wiki/De_Boor%27s_algorithm#The_algorithm
(defn- deboor
  "Evaluate the B-spline at `x` for degree `p`, `control-points`, and `knots`.
  `f+` is a function that adds two points together;
  `f*` is a function that multiplies a scalar and a point."
  [f+ f* control-points knots p x k]
  (let [tₖ₊₁ (drop (inc k) knots)]
    (loop [d (->> control-points (drop (- k p)) (take (inc p)) seq)
           tᵢ (drop (- k p -1) knots)]
      (if (next d)
        (recur (step f+ f* d x tₖ₊₁ tᵢ) (rest tᵢ))
        (first d)))))

(defn- find-index-of-knot-interval
  "Return the index of the knot interval containing `x`."
  [knots x]
  (loop [i 0
         j nil
         knots (seq knots)]
    (let [l (first knots)
          r (second knots)]
      (cond
        (nil? l) nil
        (nil? r) j
        (< x l) j
        (< x r) i
        (> x r) (recur (inc i) nil (next knots))
        :else (recur (inc i) (or j i) (next knots))))))

(comment
  (assert (nil? (find-index-of-knot-interval nil nil)))
  (assert (nil? (find-index-of-knot-interval [2] 2)))
  (assert (== (find-index-of-knot-interval [1 1] 1) 0))
  (assert (== (find-index-of-knot-interval [1 2] 1) 0))
  (assert (== (find-index-of-knot-interval [1 2] 2) 0))
  (assert (== (find-index-of-knot-interval [0 0 0 0 1 1 1 1] 1) 3))
  (assert (== (find-index-of-knot-interval [0 1 2 3 4 5 6 7] 7) 6))
  (assert (== (find-index-of-knot-interval [0 1 3 4 4 4 5 6] 4) 5))
)

(defn b-spline
  "Return the parametric B-spline for `control-points`, and `knots`.
  `f+` is a function that adds two points together;
  `f*` is a function that multiplies a scalar and a point."
  [f+ f* control-points knots]
  (assert (pos? (count control-points)))
  (let [p (- (count knots) (count control-points) 1)] ; Degree
    (assert (>= p 0))
    (fn [x]
      (when-let [k (find-index-of-knot-interval knots x)]
        (if (<= p k (- (count knots) p 1))
          (deboor f+ f* control-points knots p x k)
          (naive f+ f* control-points knots p x))))))
